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Abstract 

We calculate the one-loop electroweak corrections to e + e~ — > W + W~Z and 
e + e~ — > ZZZ and analyse their impacts on both the total cross section and some 
key distributions. These processes are important for the measurements of the quartic 
couplings of the massive gauge bosons which can be a window on the mechanism of 
spontaneous symmetry breaking. We find that even after subtracting the leading 
QED corrections, the electroweak corrections can still be large especially as the 
energy increases. We compare and implement different methods of dealing with 
potential instabilities in the routines pertaining to the loop integrals. For the real 
corrections we apply a dipole subtraction formalism and compare it to a phase-space 
slicing method. 



1 Introduction 



The LHC has just started running again and seems now to be on course for what it has 
been built for: discovery of the last remaining particle of the much successful Standard 
Model, the Higgs boson. It may well be that before this particle is uncovered we will have 
seen clear signs of New Physics that better encompasses an elementary Higgs boson. The 
conventional Minimal Supersymmetric Standard Model is the most popular example of 
such a scenario. It may however happen that the mechanism of electroweak symmetry 
breaking will remain elusive and that one has to look for subtle deviations in standard 
processes. Because of its clean environment a linear collider might be more suited for this 
purpose. 

From this perspective the study of e + e~ — > W + W~Z and e + e~ — > ZZZ may be very 
instructive and would play a role similar to e + e~ — > W + W~ at lower energies. Indeed 
it has been stressed that e + e~ — > W + W~Z and e + e" — > ZZZ are prime processes for 
probing the quartic vector boson couplings [1]. In particular deviations from the gauge 
value in the quartic W + W~ Z Z and ZZZZu couplings that are accessible in these reac- 
tions might be the residual effect of physics intimately related to electroweak symmetry 
breaking. Since these effects can be small and subtle, knowing these cross sections with 
high precision is mandatory. This calls for theoretical predictions taking into account 
loop corrections. 

Apart from the physics motivations for performing such calculations, the other reason 
is that one-loop corrections, in particular the electroweak corrections, for such 2 — > 3 
processes are a good testing ground for the various ingredients and techniques that enter 
such one-loop multi-leg corrections. Although recently NLO corrections to 2 — > 4 pro- 
cesses have set the technical frontier with a handful of processes in this category having 
been addressee!!, NLO corrections to 2 — > 3 processes are far from straightforward. Not 
only the number of diagrams differs greatly from one process to another but perhaps more 
importantly the loop structure can also differ significantly. For the processes at hand one 
has to deal with high rank tensors, rank-4, for the pentagon diagrams, compared to at 
most a rank-2 for e + e~ — > vvH [5J [6], [TJ, [8] . This might lead to much more severe numerical 
instabilities due to the appearance of higher powers of the inverse Gram determinants in 
the tensor reduction. Moreover different scales and masses may lead to sensitive issues 
related to Landau singularities in scalar integrals [H] . It is therefore important to conduct 
one-loop corrections to a variety of 2 — > 3 processes. 

Radiative corrections to e + e _ —> ZZZ have appeared recently in [10] and those to 
e + e~ — > W + W~Z in [TT] while we were preparing this paper. We have made an indepen- 
dent calculation of the electroweak corrections to e + e~ — > W + W~Z and e + e~ — > ZZZ . 
Preliminary results on e + e~ — > W + W~Z have been presented in [J_2] before those of [TT] 
were made public. We perform two independent calculations and check further through 
non-linear gauge parameter independence tests. These help also identify potential in- 
stabilities in the routines pertaining to the loop integrals. We detail how some critical 

X ZZZZ is absent at tree-level in the Standard Model. Other photonic quartic couplings may be 
probed in the processes we study but they are best studied in other reactions. Moreover the latter, from 
the point of view symmetry breaking, are less relevant and would be of higher order [2J. 

2 As far as electroweak corrections are concerned only two such calculations have been performed [3l|4]. 
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Figure 1: Representative Born diagrams for e + e~ — > ZZZ and e + e~ — > W + W~Z . Di- 
agrams (a) contribute to both processes while diagrams of type (b) contribute only to 
e + e~ — > W + W~Z . The first diagram of type (a) will be referred to as the Higgsstrahlung 
contribution. 

issues related to inverse Gram determinants have been tackled, how real corrections have 
been implemented and how checks were conducted. Our calculation is implemented in 
two independent Monte Carlo codes which can calculate total cross sections and arbitrary 
distributions. We compare our results with those in [TOT ITT] and comment also on the 
renormalisation scheme and input parameters. 

2 Calculational details 

At leading order W + W~Z and ZZZ final states are produced through the diagrams 
shown in Fig. [TJ A common contribution to the two processes is the Abelian-like t- 
channel fermionic exchange akin to the QED process e + e~ — > 37. Both processes include 
the Higgsstrahlung contribution where the splitting H* — > VV occurs. Apart from this 
contribution e + e~ — > W + W~Z can be built up from e + e~ — > W + W~ through the addi- 
tion of Z radiation from either the initial or final state and the s-channel quartic WWZZ. 
Since the precision electroweak data suggest a Higgs mass below the WW threshold, we 
restrict our study to the region Mh < 160 GeV. This means that the Higgsstrahlung 
contribution can not be resonant and therefore in our calculation no width is introduced. 

We also set the electron mass to zero whenever possible. The electron mass then 
appears only in mass singular logarithms. These arise in the virtual corrections from 
loop diagrams containing electrons and from photons radiated off electrons in the real 
corrections. 

We have performed our calculation in at least two independent ways both for the 
virtual and the real corrections leading to two independent numerical codes. A comparison 
of both codes has shown full agreement at the level of the integrated cross sections as well 
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as all the distributions that we have studied. 

The phase-space integration is done by using the Monte Carlo integrator BASES [J3J [33] 
in one code while the other code employs VEGAS [15J. 



2.1 Renormalisation 

We adopt the on-shell renormalisation scheme as detailed in JTHJ HZ]. By default, in this 
scheme the electromagnetic coupling is defined in the Thomson limit at q 2 — > 0. The 
counterterm e — > Ye = (l + 5Z e )e is related through a Ward identity to the wave function 
renormalisation constant of the photon and the wave function describing the A — > Z 
transition defined at q 2 = so that the photon propagator is defined with residue equal 
to one and no A — > Z transition remains when the photon is on shell (TJIJ 112] • In the 
conventions of [TB] this leads to 



5Z e = -5zT A + S f5Z l Jl 5Z% = i^VW 6Z£ = ~lt^(0). (1) 



Tlj. is the transverse parts of the VV self-energy. This particular definition of the charge 
at the scale q 2 = is not the most appropriate since the weak processes take place at 
scales of order Mz or higher. The running of a from q 2 = to q 2 = M§ alone amounts 
to a 6% correction. For a process of order a n , this running will amount to a correction of 
order n x 6%, thus hiding more interesting corrections. Moreover these corrections due 
to the running are sensitive to the light fermion masses through logarithms of the type 
\n(q 2 /m 2 ). In fact the effective couplings of the Z to fermions are also sensitive to isospin 
breaking effects and therefore virtual heavy top effects through Ar [18J. The combined 
effect of these two corrections is better parameterised if one uses the Fermi coupling in 
lieu of a(0). We will therefore use a variant of the scheme. This scheme absorbs a 
large universal part of 0(a) corrections into the Born contribution. The advantage of the 
scheme is that the final results are not sensitive to the light fermion masses, in particular 
the light quark masses, and some universal m 2 corrections are also absorbed. In this 
scheme we use Mz, Mw, other masses} instead of {a(0), Mz, Mw, other masses} as 
input parameters, from which the electromagnetic coupling constant is calculated as 



To avoid double counting we have to subtract the one-loop part of the universal cor- 
rection from the explicit 0(a) corrections by using the counterterm 



In the 't Hooft-Feynman gauge with the usual linear gauge fixing (Ar)!_i 00 p is given 




5Z^ = 5Z, 



e - -(Ar)i-ioop. 



(3) 
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Using ag* in tree-level calculations will therefore take care of some universal higher- 
order contributions. In the scheme the corrections are initially of order <Xq . Con- 
sidering however that the typical scale for both virtual photon exchange and real photon 
radiation is q 2 = we rescale our results so that the NLO results are of order a% a(0). 



2.2 Virtual corrections 

The virtual corrections have been evaluated using a conventional Feynman-diagram based 
approach using standard techniques in the two independent codes. The total number of 
diagrams in the 't Hooft-Feynman gauge is about 2700 including 109 pentagon diagrams 
for e + e~ — > W + W~Z and about 1800 including 64 pentagons for e + e~ — > ZZZ. This 
already shows that e + e~ — > W + W~ Z with as many as 109 pentagons is more challenging 
than e + e~ — > ZZZ. 



Code 1 

The first code uses FeynArts-3.4 [19J to generate all Feynman diagrams and amplitude 
expressions. FormCalc-6 . [201 EI] is used to simplify and generate a Fortran 77 code 
suited for the numerical evaluation of the differential cross sections. We also use SloopS 
[22| [23| 124"] an automated code that uses a few modules from FeynArts-3.4 but which 
implements the generalised non-linear gauge (NLG) [251 IT6] 



C GF = -^(fy - %taA^ - igcwPZ^W* + £ w ?-(v + 6H + i~kxz)X' 

1 IO. ry , ^ 9 ,_. , =TT ^ ,2 1 /Q 4X2 



|2 



,,. + + ^)X 3 ) 2 - ^- . (5) 

The x represents the Goldstone. We take the 't Hooft-Feynman gauge with ^ w = £z — 
£a = 1 so that no "longitudinal" term in the gauge propagators contributes. Not only does 
this make the expressions much simpler and avoids unnecessary large cancellations, but it 
also avoids the need for higher tensor structures in the loop integrals. The use of the five 
parameters, a, (3, 5, k, e is not redundant as often these parameters check complementary 
sets of diagrams. At a few random points in phase space we exploit these parameters 
to perform powerful tests on the generated matrix elements. This test can reveal for 
example numerical instabilities that are due to the reduction algorithm in some points in 
phase space. While at a regular point in phase space the non-linear gauge check attains 
a 14 digit agreement in double-precision when changing a NLG parameter from to 1, at 
non- regular points the same tests can fail. 

Five-point one-loop integrals (up to rank 4) are reduced to four-point integrals by 
using the reduction method of Denner and Dittmaier [26]. By default, four-point and 
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three-point tensor integrals are reduced to scalar integrals by using the Passarino-Veltman 
reduction algorithm [27J. With the latter we have observed serious problems of numerical 
instability related to four-point tensor integrals. This occurs when the Gram determinants 
associated to these tensor integrals, defined by det(Gs) = det(2piPj), become sufficiently 
small. We have solved this problem in two ways. In the first method, we use a simple 
extrapolation trick using the segmentation technique described in [28] when the Gram 
determinant is small enough. The condition implemented in our code is 



where p 2 max is the maximum external mass of a box diagram. In this limit the iV-point 
function of rank M is written as a combination of (iV — 1) -point functions of rank M. This 
is done directly at the level of the loop integral in momentum space before introducing any 
Feynman parameters. This implementation requires that one supplements the standard 
libraries with the reduction of the tensors of rank M = N + 1, for certain iV-point 
functions. Another way of tackling this problem in the first code is to calculate all the 
loop integrals in quadruple precision. For this study we have performed this everywhere 
in phase space and not just for the points that satisfy Eq. (jSJ). The numerical integration 
becomes very stable even in the case of very small Gram determinant. The price to 
pay is that the computation speed is about 6 times slower than using the segmentation 
method. We have obtained agreement of cross sections within integration errors between 
the segmentation method and using quadruple precision. 

We have also observed that the scalar one-loop four-point integral can show numer- 
ical problems and the library LoopTools [29, [301 EI] alone is not good enough for our 
calculations. While reverting to quadruple precision (everywhere in LoopTools) remedies 
the problem, in double precision we call other loop libraries to calculate scalar one-loop 
four-point integrals for some special cases where LoopTools fails. OneLOop [32] is used 
for some special cases with zero internal masses. Other specific cases that we have iden- 
tified are treated with DOC [33], a code to calculate scalar one-loop four-point integrals 
with complex/real masses. Generically, numerical instabilities in scalar one-loop integrals 
can originate from the following two sources. One is related to an endpoint singularity 
manifested as a pole very close to the boundary of the integration interval. The other is 
called a pinch singularity where there are two poles sitting very close to each other with 
a very small imaginary part. These are both called Landau singularities in Feynman loop 
integrals (e.g. see [3U |9] for a detailed discussion). At NLO, these singularities (up to 
four-point) are integrable but they may cause numerical instability. For the case of the 
pinch singularity the sign of the imaginary part of a pole (in the complex plane) can be 
wrongly calculated and hence a regular case where both poles sitting on the same side of 
the real axis can be numerically misidentified as a pinch singularity or vice versa. We have 
observed that in our calculations, in particular the process e + e~ — > W + W~Z, numerical 
problems in the scalar four-point integrals are related to both sources. An efficient way 
to cross check the results of scalar integrals is therefore to introduce a tiny positive width 
for internal masses (Tj = 10 -5 mj). In this case, the masses become complex and the 
code DOC must be used. This is indeed what we did in our calculations to obtain the 
preliminary results [12] which agree within integration errors with our final results. In 




(6) 
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e + e~ — > ZZZ an example where DOC is used is for = 300 GeV. Here the problem 
is related to the ti threshold (y'i < 2m t in this case) in the box diagram with four top 
quarks in the loop. Since all the important discriminants in DOC depend only on external 
momenta, the problem does not occur in DOC. 
Code 2 

The second code also uses FeynArts-3 . 4 for Fey nman- diagram and amplitude generation 
and FormCalc-6 . to evaluate the amplitudes. The analytical output of FormCalc-6 . 
in terms of Weyl-spinor chains and coefficients containing the tensor one-loop integrals is 
then translated to C++ code after performing further optimizations of the expressions. 

The evaluation of the one-loop tensor integrals is done by reducing them to a set of 
scalar integrals. The 5-point integrals are written in terms of 4-point functions follow- 
ing [33], which avoids leading inverse Gram determinants and the associated numerical 
instabilities. The remaining 3- and 4-point tensor integrals are recursively reduced to 
scalar integrals with the Passarino-Veltman algorithm [27]. For exceptional phase-space 
points this reduction scheme becomes numerically unstable. In this case we reevaluate 
both the scalar integrals and the reduction itself in higher precision using the QD library 
|36j . To determine when to switch to quadruple precision we use the condition numbeio 
of the Gram matrix. This is a good estimator of the number of digits lost in the solu- 
tion of the linear equation system appearing in Passarino-Veltman reduction. While this 
simple estimator is sufficient for triangle integrals in the case of 4- and 5-point integrals 
the numerical instabilities can also originate from small Gram determinants in the lower 
N-point integrals. We therefore use not only the condition number of the N-point Gram 
matrix but also the condition numbers of the (N-l)- and (N-2)-point integrals appearing 
in the tensor reduction in these cases. 

Finally, the scalar integrals are calculated using the results of [38] 139] for the finite 
and [lOj HI] for the IR singular integrals. The scalar integrals and the tensor reduction 
have been implemented as a C++ library allowing the calculation both in double and 
higher precision. Internally the library uses dimensional regularization for both UV and 
IR divergences. For this calculation the IR singularities are then translated to photon 
mass regularization [42] . 

2.3 Real corrections 

In addition to the virtual corrections we also have to consider real photon emission, i.e. 
the processes e + e~ — > W + W~ Zj and e + e~ — > ZZZ'j. The corresponding amplitudes are 
divergent in the soft and collinear limits. The soft singularities cancel against the ones 
in the virtual corrections while the collinear singularities are regularized by the physical 
electron mass. To extract the singularities from the real corrections and combine them 
with the virtual contribution we apply both the dipole subtraction scheme and a phase- 
space slicing method. 

The dipole subtraction formalism is a process independent approach and was originally 
introduced for massless QCD [43J. We use a generalization of this method to include 
photon radiation from massive fermions [44J. Since photon radiation from a massive 

3 The condition number of a symmetric matrix is denned as the ratio of the largest and smallest 
eigenvalues, see e.g. [57] . 
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fermion and from a massive charged gauge boson has the same singular structure in the 
soft limit this formalism can directly be applied to our calculation. 

In the dipole subtraction method a specially constructed function is subtracted from 
the real amplitude and then added back 

Oreai = / (dcr rea i - dcr sub ) + / da sub . (7) 

J 4 </4 

The subscript 4 refers to the 4-body final state including photon radiation. The subtrac- 
tion function has the same singular structure as the amplitude pointwise in phase space. 
The difference of the real amplitude and the subtraction function is therefore regular and 
the integration of the first term in Eq. (j7J) can be performed numerically. Introducing a 
photon mass as a regulator the subtraction function can be integrated analytically over 
the photon emission phase space up to a convolution 

/ d<7 sub = J dx^QiQjQijix) / dfXBorn + ^endpoint, (8) 

where the charges Qi are counted as outgoing and the subscript 3 refers to the 3-body 
phase space without photon radiation. While the first term in Eq. ([8]) has only mass 
singularities, the endpoint contribution contains both soft and collinear singularities and 
is given by 

^endpoint = — / do"Born ^ ] QiQj Gij. (9) 

271 J3 TT! 

The summations in Eq. flB]) and Eq. (Q run over all initial and final state charged particles. 
The explicit expressions for the and the Gij{x) can be found in [ll] . The soft and 
collinear singularities cancel in the sum of the endpoint and virtual contributions. The 
mass singularities in the final result then originate only from the first term in Eq. ([H]). 

The idea of phase-space slicing is to split the phase space of the photon emission con- 
tribution into a soft, a collinear and a remaining finite part. In the soft and collinear 
regions the real amplitude approximately factorizes into universal soft and collinear func- 
tions and the Born amplitude. In addition the phase space splits into the leading order 
phase space and a soft or collinear part. The phase-space integration over the photon 
degrees of freedom can then be performed analytically resulting in infrared and mass sin- 
gular contributions. In the remaining part of phase space the amplitude is regular and 
the integration can be performed using numerical integration. 

We have implemented a two-cutoff phase-space slicing method closely following (13 
145] . The differential cross section for the real contribution is decomposed as follows 

d(T r eai = da soft (8 s ) + d(T hard (8 s ) , 

da hard (5 s ) = da co u(8 s ,8 c ) + d<7f in (8 s ,8 c ) (10) 

using the soft and collinear cutoffs S s and 8 C . The soft region is defined by E y < 5 s y/s/2 = 
AE. The collinear but non-soft region is given by {E 1 > AE, 1 — cos# 7 / < 5 C } where 9 1 f 
is the polar angle of the photon with respect to the e ± direction in the cm. frame. 
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Figure 2: Dependence of cr^J W Z7 on i/ie so/t cutoff 5 S in phase-space slicing with 
fixed 8 C = 7 ■ I0~ 4 . On/y £/ie non- singular part is shown, i.e. the IR singular ln(m^) terms 
are set to zero. The result using dipole subtraction is shown for comparison with the error 
given by the width of the band. 
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Since the approximations used in the soft and collinear regions introduce errors of 
0(S S , S c ) the cutoffs have to be chosen sufficiently small. As can be seen from Fig. [2] the 
approximation errors are below the integration errors for 5 S < 1 0~ 3 . For smaller values of 
the soft cutoff the errors grow larger due to cancellations between the soft and the hard 
contributions which both diverge as log5 s . We have similarly verified the stability with 
respect to the variation of the collinear cutoff 5 C . Fig. [2] also shows agreement between 
the slicing and dipole subtraction results, although the errors are typically a factor fO 
smaller when using the subtraction formalism. 

2.4 Defining the weak corrections 

It is well known that the collinear QED correction related to initial state radiation (1SR) 
in e + e~ processes is large. The effect due to ISR can be treated along a structure function 
approach which resums the effects of higher orders, see for example jlT] . Likewise for the 
linear collider the effect of beamstrahlung [18] needs to be convoluted over the genuine 
weak corrections. Therefore in a NLO computation such as ours it is best to subtract 
the ISR corrections in order to sum their effect to all orders or, put another way, once 
the weak correction has been defined to convolute its result within a structure function 
approach. Deconvolution is also possible from the experimental data to arrive at the non- 
ISR result as is done at LEP, for example jl?]. The weak corrections encapsulate more 
interesting features. In order to see the effect of the weak corrections, one should separate 
the large QED corrections from the full NLO result. It means that we can define the weak 
correction as an infrared and collinear finite quantity. In our work we will subtract all of 
the QED corrections, not only the initial but, for WWZ, also the final and the interference 
QED corrections. The definition we adopt in this paper is based on the dipole subtraction 
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a V+S = — ( ( L e - !) + ~ L e ~ 1 ) ^Born , L e = h\{s/m 2 e ) . (12) 



formalism. In this approach, the sum of the virtual and endpoint contributions satisfies 
the above conditions^ and can be chosen as a definition for the weak correction 

'-'"weak "virt ^midpoint • (11) 

For the numerical results shown in Section [3] and Section HI we will make use of this 
definition. 

If one uses the phase-space slicing approach, a definition for the weak correction can 
be found as well. In the neutral process e + e~ — > ZZZ the QED corrections are confined 
to the initial state. The universal leading QED part of cr virt + cr so f t can be extracted and 
is given b}@ 

— (L e -l)\n5 s + 

7T V 4 

Subtracting this from the sum of virtual and soft contributions we can define the weak 
corrections in phase-space slicing as 

ZZZ _ , _ QED /in\ 

a weak, slicing — °virt ~T~ CT so f t &V+S- 

This procedure will lead to the same result for the weak correction as obtained by simply 
taking the sum of the virtual and endpoint parts. 

For e + e~ — > W + W~ Z there is also final state radiation and its interference with the 
initial state radiation. Diagrammatically there is no unambiguous way to subtract this 
contribution in a gauge invariant way. Nonetheless after subtracting the initial state 
radiation in Eq. ( 1T2|) . the logarithms of infrared origin can be easily isolated and the weak 
correction defined by 

ffvirt + crsoft ~ cr$l D s = cr^^ sUcing + b\n5 s . (14) 

The coefficient b of In S s can be extracted by choosing two different values of 8 S (which 
should be sufficiently small). We have compared the weak correction obtained in this 
way with the one calculated by taking the sum of the virtual and endpoint parts. The 
results for the case of y/s = 500 GeV and M H = 120 GeV are: 5 d j e ^ = -7.014(5)%, 
$weak 9 = — 6.73(1)%. This is for the process e + e" — > W + W~ Z and other input parameters 
as specified in Section 12.51 

2.5 Input parameters 

We use the following set of input parameters (HI \5U\ , 

G ^ = 1.16637 x 10~ 5 GeV" 2 , a(0) = 1/137.035999679, 
M w = 80.398 GeV, M z = 91.1876 GeV, 
m e = 0.510998910 MeV, m M = 105.658367 MeV, m T = 1776.84 MeV, (15) 
m u = 66 MeV, m c = 1.2 GeV, m t = 173.1 GeV, 
m d = 66 MeV, m s = 150 MeV, m b = 4.3 GeV. 



A similar discussion was given in [46] . 

5 This definition differs by the sub-sub-leading term 2a/ir x 7r 2 /6 to what is usually taken for the 
universal initial state radiation. Adding this term would give the same result as the one based on a 
diagrammatic approach as done in [TU] for ZZZ. 
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Figure 3: Left: the total cross section for e + e~ — > ZZZ as a function of \fs for the Born, 
full 0(a) and genuine weak correction for Mh = 120 GeV. Right: the corresponding 
relative percentage corrections o'nlo/o'lo — 1- 



The masses of the light quarkgj, i. e. all but the top mass, are effective parameters adjusted 
to reproduce the hadronic contribution to the photonic vacuum polarization of [52] with 
oT 1 (M\) = 128.907. As discussed in Subsection 12.11 we use a variant of the scheme 
with a G(i at leading order leading to NLO corrections that are of 0(a^ a(0)). Using a G(i 
as coupling we calculate Ar = 3.0792 x 10~ 2 for M H = 120 GeV and Ar = 3.1577 x 10~ 2 
for Mh = 150 GeV. The Cabibbo-Kobayashi-Maskawa matrix is set to be diagonal. For 
the calculation we neglect the electron Yukawa coupling proportional to the electron mass, 
as mentioned at the beginning of Section [2J For both processes we apply no cuts at the 
level of the and Z, since these will decay. 



3 e+e" ZZZ 

As shown in Fig. [3] the tree-level cross section rises sharply once the threshold for pro- 
duction opens, reaches a peak of about 1.1 fb around a centre-of-mass energy of 600 GeV 
before very slowly decreasing with a value of about 0.9 fb at 1 TeV. Exact results are 
displayed in Table [TJ The Higgsstrahlung contribution to the total cross section is about 
10% at v/i = 600 GeV and = 1 TeV. 

The full NLO corrections are quite large and negative around threshold, —35%, de- 
creasing sharply to stabilise at a plateau around y/s = 600 GeV with —16% correction. 
The sharp rise and negative corrections at low energies are easily understood. They are 
essentially due to initial state radiation (ISR) and the behaviour of the tree-level cross 
section. The photon radiation reduces the effective centre-of-mass energy and therefore 
explains what is observed in the figure. On the other hand the genuine weak corrections, in 
the scheme, are relatively small at threshold, —7%. The magnitude of the corrections 
however increases steadily reaching a value as large as —18% at \fs = 1 TeV. These large 

6 which are the same as those used in [5T] . 
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Table 1: Cross section for e + e~ — > ZZZ at tree-level, including only the weak corrections 
and at full next-to-leading order for Mn = 120 GeV. Also shown are the relative weak 
and full NLO corrections. 



^[TeV] 


OBorn [ fb] 


O"xoeofe[fb] 




&weak [%] 


Sfull[%] 


0.3 


0.184524(4) 


0.17132(1) 


0.11916(2) 


-7.156(7) 


-35.420(8) 


0.4 


0.92437(3) 


0.83637(9) 


0.73502(9) 


-9.520(9) 


-20.484(9) 


0.5 


1.12353(1) 


0.99338(9) 


0.92820(9) 


-11.584(7) 


-17.386(8) 


0.6 


1.14203(6) 


0.9900(2) 


0.9546(2) 


-13.31(1) 


-16.41(1) 


0.8 


1.04796(8) 


0.8801(2) 


0.8786(2) 


-16.02(2) 


-16.16(2) 


1.0 


0.92962(10) 


0.7609(2) 


0.7759(2) 


-18.15(2) 


-16.54(2) 



negative corrections are typical of the electroweak Sudakov logarithms — log 2 (s/M^ / ). In 
the usual a(0) on-shell scheme this important effect would be blurred and weakened unless 
one reaches much higher energies. 

We have also studied distributions in some key kinematic variables and how they are 
affected by radiative corrections. Fig. H] shows the invariant mass Mzz, the transverse 
momentum and the rapidity yz for y/s = 500 GeV and Mh = 120 GeV. The important 
message is that the genuine weak corrections are almost just an overall rescaling of the 
leading-order distributions, in particular in the bulk of the region of phase space where 
the yield is largest. The full NLO corrections on the other hand show more structure 
with very large corrections at the edges of phase space where the cross section is smallest, 
for example the full NLO correction for Mzz > 350 GeV drops below —50%. This again 
is essentially due to hard photon radiation. This shows that the effects of New Physics 
could be discovered in a less ambiguous way after subtracting the QED corrections. 

4 e+e" -> W + W~Z 

Compared to e + e~ — > ZZZ, the cross section for e + e~ — > W + W~ Z is almost 2 orders of 
magnitudes larger for the same centre-of-mass energy. For example at 500 GeV it is about 
40 fb at tree level, compared to 1 fb for the e + e~ — > ZZZ cross section. For an anticipated 
luminosity of lab -1 , this means that the cross section should be known at the per- mil level. 
In absolute terms the Higgsstrahlung contribution is a factor 2 (due to SU (2)) larger than 
in e + e~ — > ZZZ, however its contribution to the total e + e~ — > W + W~Z cross section is 
much less than 1% for M H = 120 GeV. 

The behaviour of the e + e~ — > W + W~Z cross section as a function of energy resembles 
that of e + e~ — > ZZZ. It rises sharply once the threshold for production opens, reaches 
a peak before very slowly decreasing as shown in Fig. |5] (exact results are also displayed 
in Table [2]) . However as already discussed the value of the peak is much larger, ~ 50 fb 
at NLO, moreover the peak is reached around \fs = 1 TeV, much higher than in ZZZ. 
This explains the bulk of the NLO corrections at lower energies which are dominated by 
the QED correction, large and negative around threshold and smaller at higher energies. 
As the energy increases the weak corrections get larger reaching about —18% at \fs = 



11 




e + e"^ZZZ 
\fs=500GeV 
M H =120GeV 



do NLO /do Born -1[%] 
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Figure 4: From top to bottom: distributions for the ZZ invariant mass, the rapidity 
of the Z and the transverse momentum of the Z for e + e~ — > ZZZ at ^ = 500 GeV 
and Mh = 120 GeV. The panels on the left show the tree-level, the full NLO and the 
weak correction. The panels on the right show the corresponding relative (to the tree- 
level) percentage corrections. The distributions are obtained by entering for each event 
the corresponding observable, say , of each Z and then normalising by a factor 1/3. 
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Figure 5: Left: the total cross section for e + e~ — > W + W~Z as a function of -y/i for the 
Born, full 0(a) and genuine weak correction for M# = 120 GeV. Right: the corresponding 
relative percentage corrections o~nlo/°'lo — 1- 



Table 2: Cross section for e + e~ — > W + W~Z at tree-level, including only the weak cor- 
rections and at full next-to-leading order for M H = 120 GeV. Also shown are the relative 
weak and full NLO corrections. 



^[TeV] 


CBorn [ft>] 




o-full[M 




5 f uii[%] 


0.3 


3.27055(4) 


3.1888(3) 


2.3880(3) 


-2.500(8) 


-26.986(9) 


0.5 


39.7557(9) 


36.967(2) 


33.476(2) 


-7.014(5) 


-15.795(5) 


0.7 


55.358(3) 


49.878(6) 


47.409(6) 


-9.899(10) 


-14.359(10) 


0.9 


59.121(4) 


51.881(8) 


50.678(8) 


-12.25(1) 


-14.28(1) 


1.0 


59.061(4) 


51.206(9) 


50.541(9) 


-13.30(1) 


-14.43(1) 


1.2 


57.202(5) 


48.49(1) 


48.69(1) 


-15.24(2) 


-14.88(2) 


1.5 


52.740(5) 


43.34(1) 


44.43(1) 


-17.82(2) 


-15.76(2) 



1.5 TeV. Again a large part of this correction seems to be of the Sudakov type. 

More interesting than in the case of e + e _ — > ZZZ are the distributions in some 
key variables like the invariant WW mass, the and the rapidity of the WW system. 
First, due to photon radiation, in the full NLO corrections some large corrections do 
show up at the edges of phase space, see Fig. [61 However when the QED corrections are 
subtracted, the weak corrections cannot be parameterised by an overall scale factor, for 
all the distributions that we have studied. Other distributions not shown here can be 
found in |53j . 
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Figure 6: From top to bottom: distributions for the WW invariant mass, the rapidity 
of the WW system and the transverse momentum of the Z for e + e~ — » W + W~ Z at 
y/s = 500 GeV and M H = 120 GeV. The panels on the left show the tree-level, the full 
NLO and the weak correction. The panels on the right show the corresponding relative (to 
the tree-level) percentage corrections. 
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5 Comparison with other calculations 



The electroweak corrections to e + e~ — > ZZZ have also been calculated in [10]. They 
use the a(0) scheme and slightly different input parameters. We have performed a tuned 
comparison by adapting to the input parameters of [TO] and switching to the a(0) scheme. 
We find full agreement within the quoted statistical errors for o^q and A<j toi shown in 
Table 2 of [TU] . A comparison of the results in Fig. 4 of [10] is shown in Table [^|. The 
results for the Born cross section agree to within 0.01% while the NLO results agree to 
at least 0.1%. 

Table 3: Born cross section and relative corrections for e + e~ — > ZZZ using the input 
parameter scheme of [TUf . 







M H = 


120 GeV 


M H = : 


L50 GeV 


v^[GeV] 




CBom [fb] 


5 fun [%] 


&Born [fb] 


Sfull [%} 


350 


Ref. [TO] 


0.58696 


-15.79 


0.68422 


-13.91 




This work 


0.586955(2) 


-15.850(1) 


0.684209(2) 


-13.970(1) 


370 


Ref. [10] 


0.70531 


-13.79 


0.80821 


-12.00 




This work 


0.705303(2) 


-13.822(1) 


0.808196(3) 


-11.986(1) 


400 


Ref. [TO] 


0.83409 


-11.75 


0.9375 


-9.98 




This work 


0.834083(4) 


-11.765(2) 


0.937484(4) 


-9.973(1) 


450 


Ref. [TO] 


0.95792 


-9.79 


1.05294 


-8.06 




This work 


0.957904(5) 


-9.763(3) 


1.052917(5) 


-8.044(2) 


500 


Ref. [TO] 


1.01384 


-8.70 


1.09754 


-7.09 




This work 


1.013806(6) 


-8.682(4) 


1.097440(7) 


-7.064(4) 


600 


Ref. [TO] 


1.03052 


-7.77 


1.09370 


-6.36 




This work 


1.030489(9) 


-7.714(6) 


1.093668(9) 


-6.289(6) 


700 


Ref. [TO] 


0.99611 


-7.47 


1.04437 


-6.20 




This work 


0.99607(1) 


-7.438(9) 


1.04437(1) 


-6.164(9) 


800 


Ref. [TO] 


0.94567 


-7.50 


0.98647 


-6.61 




This work 


0.94563(1) 


-7.46(1) 


0.98343(1) 


-6.30(1) 


900 


Ref. [TO] 


0.89168 


-7.71 


0.92196 


-6.65 




This work 


0.89164(1) 


-7.62(1) 


0.92191(1) 


-6.55(1) 


1000 


Ref. [TO] 


0.83892 


-7.94 


0.86366 


-6.89 




This work 


0.83887(2) 


-7.86(2) 


0.86362(2) 


-6.86(2) 



We have also made a comparison for the process e + e~ — > W + W~ Z with the results 
of [IT] . In addition to different input parameters, [UJ uses an unusual scheme/input 
parameter for the electromagnetic coupling constant. One can read from [UJ that their 
renormalisation condition for the electric charge is the on-shell definition at q 2 = 0, see 
Eq. (P), yet the value of cnjjg(Mz) is used as a value for the coupling already at tree 
level and no shift of the electric charge counterterm at NLO is made in order to avoid 
double counting. This is just like taking a numerical value of «m^(M|) for a(0) in the 
usual on-shell scheme. For the sake of comparison we have used the same scheme in our 

7 The data for Fig. 4 of [TU] have been kindly provided by Prof. Ma. 
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calculation and show the comparison with Table 2 of [11] in Table HI The agreement for 
the Born result is within integration errors as are the NLO corrections for Mjj = 120 GeV. 
However for Mh = 150 GeV we find agreement only for yfs = 0.5 TeV. Close to threshold 
and especially at high energy the results differ by up to 5 times the integration errors. 

Table 4: Born cross section and NLO correction for e + e~ — > W + W~Z using the input 
parameter scheme of FIT] /. 







M H = 


120 GeV 


M H = 


150 GeV 


v^[TeV] 




&Bom [fb] 


Aa NLO [fb] 


C Born [ fb] 


Aa NLO [fb] 


0.3 


Ref. [TT] 


3.6216(2) 


-0.683(2) 


3.8856(2) 


-0.694(2) 




This work 


3.62165(5) 


-0.6901(3) 


3.88558(5) 


-0.7010(3) 


0.5 


Ref. [Ti] 


44.026(5) 


-3.03(6) 


44.303(5) 


-2.89(6) 




This work 


44.0235(10) 


-3.107(3) 


44.301(1) 


-2.949(3) 


0.8 


Ref. [TI] 


64.35(1) 


-3.48(7) 


64.50(1) 


-3.57(9) 




This work 


64.345(4) 


-3.466(8) 


64.488(4) 


-3.250(8) 


1.0 


Ref. p] 


65.42(1) 


-3.74(9) 


65.51(1) 


-3.90(9) 




This work 


65.401(5) 


-3.650(9) 


65.499(5) 


-3.440(10) 



Note Added 

After our paper was made public, [H] was updated. The on-shell a(0) scheme is now 
implemented properly. Moreover they have improved the numerical stability of the phase- 
space integration at high energies and the results for WWZ at y/s = 800 GeV and 1 TeV 
have substantially changed. We now find agreement for these energies while the small 
discrepancy near threshold (^/i = 300 GeV) remains. 

6 Conclusions 

We have performed a calculation of the full next-to-leading order correction to the pro- 
cesses e + e~ — > W + W~ Z and e + e~ — > ZZZ in the energy range of the international linear 
collider and for Higgs masses below the WW threshold. These processes would be the 
successor of e + e~ —> W + W~ in that they would measure the quartic couplings WW ZZ 
and ZZZZ which could retain residual effects of the physics of electroweak symmetry 
breaking. With this in mind we have subtracted the QED corrections and studied the 
genuine weak corrections in the scheme. We find that the weak corrections can be 
large. For example, for a centre-of-mass energy of 1 TeV these corrections reach —13% for 
WWZ and —18% for ZZZ and grow larger for higher energies. At lower energies around 
the production threshold the cross sections are small and the weak corrections are modest. 
However, in this energy range the QED corrections are largest due to the rapid rise of the 
cross section. We have also studied the effects of the genuine weak radiative corrections 
on various distributions. While for the ZZZ channel the effect might be described by an 
overall rescaling over most of the range of the kinematic variable under consideration, the 
corrections in the WWZ channel show more structure, pointing once again to the need 
for radiative corrections when looking for New Physics effects. 
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The calculations involved in our studies are also a contribution to the field of one- 
loop corrections for multi-leg processes and the techniques that one requires to control 
and further develop to make these calculations as automatic as possible. Since these 
processes involve three-vector boson production they are much more complex than say 
single Higgs production not only from just a counting of the one- loop diagrams involved 
but also the fact that the loop integrals are more challenging, in particular for e + e~ — > 
W + W~ Z. We have shown that using higher precision arithmetic or a combination of 
different loop integral libraries can be very efficient to overcome numerical instabilities. 
These techniques are instrumental for a very precise computation, especially for the WWZ 
channel where the foreseen luminosity calls for a better than per-mil accuracy for the 
theoretical prediction. We have shown that our results agree well with those of a previous 
calculation of e + e~ — > ZZZ in [JO] but not as well, especially for some Higgs masses and 
centre-of-mass energies, for the WWZ channel in [TTj . 
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